Introduction

With this set of experiments, I aimed to build SLiM models that tracked the evolution of 2 traits. The SLiM script uses one mutation type for all QTL mutations (m3) which have their effects pulled from a multivariate normal distribution, with mean = 0 for both traits, and a variance covariance matrix (sigma) with variances of both traits being 1.0 and the covariance between the traits being varied between treatments (from 0.125 to 0.5). This covariance represents the pleiotropic effects of the mutations on the two traits, with a covariance value independent from the effects of linkage.

I plotted the distribution of population means for both traits with the various treatments of pleiotropic effects, as well as the previously-included parameters of deleterious mutations (delmu) and recombination/linkage (r).

As well as this, I did it for population variances for both traits, and the covariances.

The SLiM Script

The most important parts are at the end of initialize(), in which we set up parameters for the distribution we’re pulling phenotype effects from. Also, mutation() callback where we draw those effects when mutations arise and store them as a tagged value to be pulled in the final output (in late()).

// set up a simple neutral simulation, mu = mutation rate, Ne = eff pop size, varmu = mutational variance-
initialize() {
    setCfgParam("seed", asInteger(runif(1, 1, 2^62 - 1)));
    setSeed(seed);          //for command line, multiple runs
    setCfgParam("mu", 3.24e-9);
    setCfgParam("Ne", 500);
    setCfgParam("del_mean", -0.03);
    setCfgParam("del_shape", 0.2);
    setCfgParam("delmu", 0.1);
    setCfgParam("recom", 1e-8);
    setCfgParam("r", 1e-8);
    
    initializeMutationRate(mu);
    
    initializeMutationType("m1", 0.5, "f", 0.0);
    
    //m2 mutation type: deleterious/background selection
    initializeMutationType("m2", 0.5, "g", del_mean, del_shape);  // del_mean and del_shape are defined in the batch run script, testonemultipleruns.sh
    
    //m3 mutation type: QTL for all traits, DPE defined separately (distribution of phenotype effects)
    initializeMutationType("m3", 0.5, "f", 0.0);
    m3.convertToSubstitution = F;
    m3.color = "green";
    
    // g1 genomic element type: uses m1 and m2 for all mutations
    initializeGenomicElementType("g1", c(m1, m2), c(1, delmu));
    
    //g2 genomic element type: uses m1, m2, and m3 for all mutations
    initializeGenomicElementType("g2", c(m1, m2, m3), c(1, delmu, 1));
    
    
    // Chromosome set up with 10 loci (g2)
    initializeGenomicElement(g1, 0, 19999);
    initializeGenomicElement(g2, 20000, 29999);
    initializeGenomicElement(g1, 30000, 39999);
    initializeGenomicElement(g2, 40000, 49999);
    initializeGenomicElement(g1, 50000, 59999);
    initializeGenomicElement(g2, 60000, 69999);
    initializeGenomicElement(g1, 70000, 79999);
    initializeGenomicElement(g2, 80000, 89999);
    initializeGenomicElement(g1, 90000, 99999);
    initializeGenomicElement(g2, 100000, 109999);
    initializeGenomicElement(g1, 110000, 119999);
    initializeGenomicElement(g2, 120000, 129999);
    initializeGenomicElement(g1, 130000, 139999);
    initializeGenomicElement(g2, 140000, 149999);
    initializeGenomicElement(g1, 150000, 159999);
    initializeGenomicElement(g2, 160000, 169999);
    initializeGenomicElement(g1, 170000, 179999);
    initializeGenomicElement(g2, 180000, 189999);
    initializeGenomicElement(g1, 190000, 199999);
    initializeGenomicElement(g2, 200000, 209999);
    initializeGenomicElement(g1, 210000, 219999);
    
    rates = c(r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r, recom, r);
    ends =  c(19999, 20000, 29999, 30000, 39999, 40000, 49999, 50000, 59999, 60000, 69999, 70000, 79999, 80000, 89999, 90000, 99999, 100000, 109999, 110000, 119999, 120000, 129999, 130000, 139999, 140000, 149999, 150000, 159999, 160000, 169999, 170000, 179999, 180000, 189999, 190000, 199999, 200000, 209999, 210000, 219999);
    initializeRecombinationRate(rates, ends);
    
    // QTL constants
    
    setCfgParam("QTL_mu", c(0, 0));
    setCfgParam("QTL_cov", 0.0);
    setCfgParam("QTL_sigma", matrix(c(1, QTL_cov, QTL_cov, 1), nrow=2));
    
    catn("\nQTL DFE means: ");
    print(QTL_mu);
    catn("\nQTL DFE variance-covariance matrix (M matrix): ");
    print(QTL_sigma);
}

// This M matrix defines the amount of pleiotropy - covariation at the end includes effects
// of linkage, will likely not be this value exactly, as this is independent of other effects

// Function to make parameters work in SLiMGUI as well as command line
function (void) setCfgParam (s$ name, ifls value)
{
    if (!exists(name))
        defineConstant(name, value);
}


// create a population of Ne individuals
1 {
    sim.addSubpop("p1", Ne);
}

mutation(m3) {
    // Draw mutational effects for the m3 mutation
    effects = rmvnorm(1, QTL_mu, QTL_sigma);
    mut.setValue("e0", effects[0]);
    mut.setValue("e1", effects[1]);
    
    return T;

}


1:20000 late() {
    
    if (sim.generation % 100 != 0)
        return;
    
    for (ind in sim.subpopulations.individuals)
    {
        muts = ind.genomes.mutationsOfType(m3);
        mutscount = size(muts);
        
        phenotype0 = mutscount ? sum(muts.getValue("e0")) else 0.0;
        phenotype1 = mutscount ? sum(muts.getValue("e1")) else 0.0;
        ind.setValue("phenotype0", phenotype0);
        ind.setValue("phenotype1", phenotype1);
    }
    
    inds = sim.subpopulations.individuals;
    pheno0mean = mean(inds.getValue("phenotype0"));
    pheno1mean = mean(inds.getValue("phenotype1"));
    pheno0var = var(inds.getValue("phenotype0"));
    pheno1var = var(inds.getValue("phenotype1"));
    phenocov = cov(inds.getValue("phenotype0"), inds.getValue("phenotype1"));
    
    meanline = paste(c(sim.generation, ", ", pheno0mean, ", ", pheno0var, ", ", pheno1mean, ", ", pheno1var, ", ", phenocov, ", ", getSeed(), ", ", recom, ", ", delmu, ", ", QTL_cov, ""));
    file = paste(meanline, "");
    
// SLiMGUI output   
//  catn(paste(c("Generation: ", sim.generation, " Phenotype 0 Mean: ", pheno0mean, " Phenotype 1 Mean: ", pheno1mean, " covariance: ", phenocov)));

// File output
    writeFile("out_Neu2T10L.csv", file, append = T);

}

20001 late() {
    sim.simulationFinished();
}

We output means of both phenotypes, their variances, and covariances

Mean densities

Using plotly to create 3D figures of each treatment combination’s density plot, showing their multivariate distribution and skew at generation 20000. Density comes from the MASS package, using kde2d (kernel distribution estimation for 2 dimensions). Need to find another function to do more than bivariate distributions.

slim_out <- read.csv("out_Neu2T10L.csv") 

names(slim_out) <- c("gen", "mean0", "var0", "mean1", "var1", "cov", "pear", "seed", "r", "delmu", "pleio")

sbstfinal <- subset(slim_out, gen = 20000)


k <- 1
plotslist <- htmltools::tagList() 

for (n in 1:length(unique(sbstfinal$pleio))) {
  for (i in 1:length(unique(sbstfinal$delmu))) {
    for (j in 1:length(unique(sbstfinal$r))) {
        
      den3d_temp <- kde2d(sbstfinal$mean0[sbstfinal$pleio == unique(sbstfinal$pleio)[n] & 
                                                      sbstfinal$delmu == unique(sbstfinal$delmu)[i] &
                                                      sbstfinal$r == unique(sbstfinal$r)[j]], 
                          sbstfinal$mean1[sbstfinal$pleio == unique(sbstfinal$pleio)[n] & 
                                          sbstfinal$delmu == unique(sbstfinal$delmu)[i] &
                                          sbstfinal$r == unique(sbstfinal$r)[j]], n = 100)
      
      
      plotmeans_temp <- plot_ly(x = den3d_temp$x, y = den3d_temp$y, z = den3d_temp$z) %>% add_surface() 
      plotmeans_temp <- plotmeans_temp %>% layout(
        title = paste0("Density of phenotypes where pleiotropy covariation = ", unique(sbstfinal$pleio)[n], "\n Proportion of deleterious mutations = ", unique(sbstfinal$delmu)[i], "\n and recombination between regions = ", unique(sbstfinal$r)[j]),
        scene = list(
          xaxis = list(title = "Phenotype 0 mean"),
          yaxis = list(title = "Phenotype 1 mean"),
          zaxis = list(title = "Density")
        ))

      plotslist[[k]] <- plotmeans_temp
      k <- k + 1
      rm(den3d_temp, plotmeans_temp)
      }
    }
  }

Results

Pleiotropy 0.125: * Deleterious mutations 0, recombination 0: One large peak in density at pheno 0 mean (p0m) = -1.96 and pheno 1 mean (p1m) = 1.52. Not much variation around this, this spike has a density ~ 0.4. Suggests that when there are no deleterious mutations to randomly remove certain alleles and complete linkage so that regions are essentially one big block that moves together, the mean phenotype values don’t stray too far between seeds, suggesting limited random effects that affect the mean - because all the mutations are all definitely going to affect the mean, their individual effect is smaller, so differences between individuals (or runs in this case) are miniscule. + when recombination is 1e-8, similar result but even greater peak density at ~0.95 - could be because with limited recombination between regions, not all mutations will impact runs in the same way - by chance some alleles in separate regions could be separated through recombination, reducing the relative effects of some alleles on the mean, as they exist in combination with various others in some individuals but not others. However this only happens some of the time, so individually, most alleles have tiny effects on phenotype, made even smaller in the case of the ones which have been recombined. + When recombination is 0.5 (completely unlinked regions), the curve fattens a bit, with the peak density being only 0.08. This is because now that each region is completely separated, between individuals it is highly likely that individuals will have completely different configurations of gene regions, so the variance in phenotype should be higher, even if each region only contributes a small change as explained above.

Pleiotropy 0.25: * Deleterious mutations 0, recombination 0: similar to pleiotropy 0.125 for the same treatment, p0m and p1m are -1 and 0.97 respectively, so clearly the effects of evolutionary forces are still far more important than those of pleiotropy. + Similar result again for 1e-8, higher density and not much variation around it. + Same for 0.5, refer to results for pleiotropy 0.125. The difference is that the means are slightly closer together, probably due to the higher pleiotropic covariance term. * Deleterious mutations 0.1, recombination 0: similar to the 0.125 results again, however there is slightly less variation around the peak, and the peak density itself is higher at ~0.4. Why? Could be just random chance - wouldn’t expect pleiotropy to affect the peaks themselves, only the covariation between the means. + Again, peak is shorter than r = 0.0 just like when pleiotropy = 0.125, although it isn’t as short (the was higher though) + Very different to the 0.125 figure - much more variation, lower peak etc.

Pleiotropy 0.5: * Deleterious mutations 0, recombination 0:

Variance through time

Here we plot variance of each mean through time, for each treatment.

## Brownian walks of variances

seedchoice <- sample(unique(slim_out$seed), 2)

for (n in 1:length(unique(slim_out$pleio))) {
  for (i in 1:length(unique(slim_out$delmu))) {
    for (j in 1:length(unique(slim_out$r))) {
      plotvar0_temp <- ggplot(subset(slim_out, pleio == unique(slim_out$pleio)[n] & 
                                    delmu == unique(slim_out$delmu)[i] &
                                    r == unique(slim_out$r)[j]), 
                   aes(x = gen, y = var0, colour = as.factor(seed))) +
              geom_line() +
              gghighlight(seed == as.factor(seedchoice[1]) | seed == as.factor(seedchoice[2])) +
              ggtitle(paste0("Walk of phenotype 0 variance over time with pleiotropic covariance = ", 
                             unique(slim_out$pleio)[n], ", deleterious mutation prevalence = ", 
                             unique(slim_out$delmu)[i], ",\nand recombination between regions = ", 
                             unique(slim_out$r)[j])) +
              theme_classic() +
              theme(legend.position = "none") +
              labs(x = "Generation", y = "Phenotype 0 variance")
      
      
      plotvar1_temp <- ggplot(subset(slim_out, pleio == unique(slim_out$pleio)[n] & 
                                       delmu == unique(slim_out$delmu)[i] &
                                       r == unique(slim_out$r)[j]), 
                              aes(x = gen, y = var1, colour = as.factor(seed))) +
        geom_line() +
        gghighlight(seed == as.factor(seedchoice[1]) | seed == as.factor(seedchoice[2])) +
        ggtitle(paste0("Walk of phenotype 1 variance over time with pleiotropic covariance = ", 
                       unique(slim_out$pleio)[n], ", deleterious mutation prevalence = ", 
                       unique(slim_out$delmu)[i], ",\nand recombination between regions = ", 
                       unique(slim_out$r)[j])) +
        theme_classic() +
        theme(legend.position = "none") +
        labs(x = "Generation", y = "Phenotype 1 variance")
      
      print(plotvar0_temp + plotvar1_temp + plot_layout(nrow = 2))
      rm(plotvar0_temp, plotvar1_temp)
    }
  }
}

Covariance through time

And here is the covariance between the traits

seedchoice <- sample(unique(slim_out$seed), 2)

for (n in 1:length(unique(slim_out$pleio))) {
  for (i in 1:length(unique(slim_out$delmu))) {
    for (j in 1:length(unique(slim_out$r))) {
      plotcov_temp <- ggplot(subset(slim_out, pleio == unique(slim_out$pleio)[n] & 
                                       delmu == unique(slim_out$delmu)[i] &
                                       r == unique(slim_out$r)[j]), 
                              aes(x = gen, y = cov, colour = as.factor(seed))) +
        geom_line() +
        gghighlight(seed == as.factor(seedchoice[1]) | seed == as.factor(seedchoice[2])) +
        ggtitle(paste0("Walk of phenotype covariance over time with pleiotropic covariance = ", 
                       unique(slim_out$pleio)[n], ",\ndeleterious mutation prevalence = ", 
                       unique(slim_out$delmu)[i], ",\nand recombination between regions = ", 
                       unique(slim_out$r)[j])) +
        theme_classic() +
        theme(legend.position = "none") +
        labs(x = "Generation", y = "Covariance between phenotypes 0 and 1")
      print(plotcov_temp)
      rm(plotcov_temp)
      
    }
  }
}

Results

As pleiotropy increases, there is an increased bias towards positive covariances, as expected. Higher Deleterious mutation rates can cause more frequent spikes in covariance as they take out hitchhiking allelic variation. Recombination doesn’t seem to do much by itself, but interacts with deleterious mutation rate to exacerbate or relieve its effects: intermediate rates